{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 磁性质数值导数 (1)：RHF 的非 GIAO 磁化率"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 创建时间：2020-08-27"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这篇文档中，我们会讨论使用 PySCF 以及其作为 libcint 的接口，计算非 GIAO 的 RHF 数值磁化率的程序。该文档大量参考 PySCF 的代码 [magnetizability/rhf.py](https://github.com/pyscf/pyscf/blob/master/pyscf/prop/magnetizability/rhf.py) 与 [nmr/rhf.py](https://github.com/pyscf/pyscf/blob/master/pyscf/prop/nmr/rhf.py)。一些公式符号参考 Atkins, Friedman [^Atkins-Friedman.Oxford.2010]。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们的讨论中所使用到的分子体系 `mol` 会是非对称的氨分子，并且取用最小基组。规范原点 (Gauge Origin) 会取在坐标原点上 `coord_orig`。其 RHF 计算放在实例 `mf`，而磁化率计算实例会放在 `mf_mag`。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import gto, scf\n",
    "from pyscf.prop import nmr, magnetizability\n",
    "import numpy as np\n",
    "\n",
    "np.set_printoptions(precision=5, linewidth=150, suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol = gto.Mole()\n",
    "mol.atom = \"\"\"\n",
    "N  0.  0.  0.\n",
    "H  0.  1.  0.2\n",
    "H  0.1 0.3 1.5\n",
    "H  0.9 0.4 -.2\n",
    "\"\"\"\n",
    "mol.basis = \"STO-3G\"\n",
    "mol.verbose = 0\n",
    "mol.build()\n",
    "coord_orig = np.zeros(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其自洽场能量为"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-55.253540514686556"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf = scf.RHF(mol).run()\n",
    "mf.e_tot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其磁化率张量 $\\xi_{ts}$ 为 (其中，$t, s \\in \\{ x, y, z \\}$ 表示三个坐标方向，需要注意这里选择了规范原点为坐标原点，若选取其它坐标则会得到非常不同的结果)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-4.94475,  0.21773, -0.08268],\n",
       "       [ 0.21773, -4.27801,  0.49885],\n",
       "       [-0.08268,  0.49885, -4.15348]])"
      ]
     },
     "execution_count": 143,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_mag = magnetizability.RHF(mf)\n",
    "mf_mag.gauge_orig = coord_orig\n",
    "mf_mag.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 基础概念"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 分子能量作为外加微扰量的函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们指出，磁化率可以看作是分子处在某一恒定外加磁场 $\\boldsymbol{\\mathscr{B}}$ 下 (作为三维矢量)，所产生的能量变化的表征：\n",
    "\n",
    "$$\n",
    "E_\\mathrm{tot} (\\boldsymbol{\\mathscr{B}}) = E_\\mathrm{tot}^{(0)} + E_\\mathrm{tot}^{(1)} \\boldsymbol{\\mathscr{B}} + E_\\mathrm{tot}^{(2)} \\boldsymbol{\\mathscr{B}}^2 + \\cdots\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以一般物理的约定俗成而言，对于外加的磁场微扰 (Atkins and Friedman, eq 13.34)\n",
    "\n",
    "$$\n",
    "E_\\mathrm{tot}^{(2)} = - \\frac{1}{2} \\boldsymbol{\\mathscr{B}}^\\dagger \\boldsymbol{\\xi} \\boldsymbol{\\mathscr{B}} = - \\frac{1}{2} \\sum_{t, s \\in \\{ x, y, z \\}} \\mathscr{B}_t \\xi_{ts} \\mathscr{B}_s\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，$\\boldsymbol{\\xi}$ 是二维对称矩阵 (或称张量，如之前代码所展示)。我们使用了 $\\boldsymbol{\\xi}$ (Atkins and Friedman, eq 13.34, termed as *magnetizability*) 而非 $\\boldsymbol{\\chi}$ (Atkins and Friedman, eq 13.3c, termed as *magnetic susceptibility*) 来表示磁化率。因此，磁化率本身可以表示为 (矩阵元的形式与向量形式)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\xi_{ts} = - \\frac{\\partial^2 E_\\mathrm{tot} (\\boldsymbol{\\mathscr{B}})}{\\partial \\mathscr{B}_t \\partial \\mathscr{B}_s},\n",
    "\\quad \\boldsymbol{\\xi} = - \\boldsymbol{\\nabla}_{\\boldsymbol{\\mathscr{B}}} \\boldsymbol{\\nabla}_{\\boldsymbol{\\mathscr{B}}}^\\dagger E_\\mathrm{tot} (\\boldsymbol{\\mathscr{B}})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 哈密顿算符作为外加微扰量的算符"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "能量可以通过波函数在哈密顿算符的变分极小值处的期望获得：\n",
    "\n",
    "$$\n",
    "E_\\mathrm{tot} (\\boldsymbol{\\mathscr{B}}) = \\langle \\Psi (\\boldsymbol{\\mathscr{B}}) | \\hat H (\\boldsymbol{\\mathscr{B}}) | \\Psi (\\boldsymbol{\\mathscr{B}}) \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，\n",
    "\n",
    "$$\n",
    "\\hat H (\\boldsymbol{\\mathscr{B}}) = \\sum_{i} \\hat h (\\boldsymbol{\\mathscr{B}}, \\boldsymbol{r}_i) + \\hat V_\\mathrm{ee} + \\hat V_\\mathrm{NN}\n",
    "$$\n",
    "\n",
    "上述算符是体系的多电子总哈密顿算符；而 $\\hat h (\\boldsymbol{\\mathscr{B}})$ 则是单电子的 Core Hamiltonian 算符；$\\hat V_\\mathrm{ee}$ 为电子互斥算符，$\\hat V_\\mathrm{NN}$ 为原子核互斥算符。需要注意，由于我们不使用 GIAO，因此 $\\hat V_\\mathrm{ee}$ 就是普通的电子互斥算符，不受外场 $\\boldsymbol{\\mathscr{B}}$ 干扰；但使用 GIAO 的情况下，可能需要额外考虑这部分贡献。\n",
    "\n",
    "$$\n",
    "\\hat h (\\boldsymbol{\\mathscr{B}}) = \\hat h {}^{(0)} + \\hat h {}^{(1)} (\\boldsymbol{\\mathscr{B}}) + \\hat h {}^{(2)} (\\boldsymbol{\\mathscr{B}})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\hat h {}^{(0)}$ 是没有外加场的算符 (这与自洽场计算过程所用到的算符相同)。其余的算符则为 (Atkins and Friedman, eq 13.26, eq 13.29)\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat h {}^{(1)} (\\boldsymbol{\\mathscr{B}}) &= \\frac{1}{2} \\boldsymbol{\\mathscr{B}} \\cdot \\boldsymbol{r} \\times \\boldsymbol{\\hat{p}} \\\\\n",
    "\\hat h {}^{(2)} (\\boldsymbol{\\mathscr{B}}) &= \\frac{1}{8} \\big( \\boldsymbol{\\mathscr{B}}^2 \\boldsymbol{r}^2 - (\\boldsymbol{\\mathscr{B}} \\cdot \\boldsymbol{r})^2 \\big)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，\n",
    "\n",
    "$$\n",
    "\\boldsymbol{r} = \\begin{pmatrix} x \\\\ y \\\\ z \\end{pmatrix}, \\quad\n",
    "\\boldsymbol{\\hat p} = \\begin{pmatrix}\n",
    "  \\displaystyle - i \\frac{\\partial}{\\partial x} \\\\\n",
    "  \\displaystyle - i \\frac{\\partial}{\\partial y} \\\\\n",
    "  \\displaystyle - i \\frac{\\partial}{\\partial z}\n",
    "\\end{pmatrix} = -i \\nabla \\boldsymbol{r}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，\n",
    "\n",
    "$$\n",
    "E_\\mathrm{tot}^{(2)} = 2 \\langle \\Psi^{(0)} (\\boldsymbol{\\mathscr{B}}) | \\sum_i \\hat h {}^{(1)} (\\boldsymbol{\\mathscr{B}}, \\boldsymbol{r}_i) | \\Psi^{(1)} (\\boldsymbol{\\mathscr{B}}) \\rangle + \\langle \\Psi^{(0)} | \\sum_i \\hat h {}^{(2)} (\\boldsymbol{\\mathscr{B}}, \\boldsymbol{r}_i) | \\Psi^{(0)} \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "前一项被称为顺磁项 (Paramagnetic)，后一项称为抗磁项 (Diamagnetic)。$\\Psi^{(0)}$ 是指未微扰的体系哈密顿算符 $\\hat H {}^{(0)}$ 的本征态，$\\Psi^{(1)} (\\boldsymbol{\\mathscr{B}})$ 则是一阶微扰的波函数；其解析的求取方法是在程序中表示为 U 矩阵，通过 CP-HF 方程求取；我们在这里不会对解析方法作说明，但了解这两项的区分是有帮助的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Core Hamiltonian 的程序实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们知道，PySCF 中，在自洽场实例中更改 Core Hamiltonian 的类方法函数 (method function) 就可以实现外场微扰下的能量计算。这在 pyxdh 偶极矩的计算 [文档](https://py-xdh.readthedocs.io/zh_CN/latest/numdiff/num_dip.html) 中有所说明。在这里我们也要做类似的工作。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 顺磁项"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在 PySCF 中，顺磁项 $\\hat h {}^{(1)} (\\boldsymbol{\\mathscr{B}}) = \\frac{1}{2} \\boldsymbol{\\mathscr{B}} \\cdot \\boldsymbol{r} \\times \\boldsymbol{\\hat p}$ 有其对应的积分 `hcore_1` ($h_{t \\mu \\nu}^{(1)}$，需要注意它不包含作为标量的 $\\mathscr{B}_t$)\n",
    "\n",
    "$$\n",
    "h_{t \\mu \\nu}^{(1)} \\cdot \\mathscr{B}_t = \\langle \\mu | \\hat h {}^{(1)} (\\mathscr{B}_t) | \\nu \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3, 8, 8), dtype('complex128'))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hcore_1 = - 0.5 * mol.intor(\"int1e_cg_irxp\") * 1j\n",
    "hcore_1.shape, hcore_1.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述的程序看起来会有些奇怪，因为这里出现了复数。我们需要分段对其作解释。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**积分字符**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们使用到了积分字符 `int1e_cg_irxp`。关于这段字符，其意义需要通过 [auto_intor.cl](https://github.com/sunqm/libcint/blob/abf6948fa17e5b4ecbd26de05bf4b1d7b2b2fe3c/scripts/auto_intor.cl#L12) 程序了解：\n",
    "\n",
    "```lisp\n",
    "  '(\"int1e_cg_irxp\"             (#C(0 1) \\| rc cross p))\n",
    "```\n",
    "\n",
    "其右侧是积分的具体形式，说明在 [README](https://github.com/sunqm/libcint/blob/master/README) 文件中，意义为\n",
    "\n",
    "$$\n",
    "\\mathtt{int1e\\_cg\\_irxp} = i \\langle \\mu | \\boldsymbol{r} \\times \\boldsymbol{\\hat p} | \\nu \\rangle\n",
    "$$\n",
    "\n",
    "其维度是 $(t, \\mu, \\nu)$，但其第一个维度是通过向量叉乘给出，因此它与 $\\boldsymbol{r}$ 或 $\\boldsymbol{p}$ 的维度不是直接相关的。如果我们令动量算符 $\\boldsymbol{\\hat l} = \\boldsymbol{r} \\times \\boldsymbol{\\hat p}$，那么可以将上述积分写为\n",
    "\n",
    "$$\n",
    "\\mathtt{int1e\\_cg\\_irxp}_{t \\mu \\nu} = i \\langle \\mu | \\hat l_t | \\nu \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**反对称性厄米性**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们应当留意到 $\\mathtt{int1e\\_cg\\_irxp}_{t \\mu \\nu}$ 是一个反对称矩阵 (即随 $\\mu, \\nu$ 交换成相反值)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(mol.intor(\"int1e_cg_irxp\"), - mol.intor(\"int1e_cg_irxp\").swapaxes(-1, -2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这是由于 $\\nabla$ 算符本身是一个反对称算符。但是需要留意到，动量算符在此基础上乘以了虚数单位 $- i$，因此，该矩阵是厄米的，即其转置后的共轭是其本身。我们所定义的 $h_{t \\mu \\nu}^{(1)}$ 就具有这样的性质："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(hcore_1, hcore_1.swapaxes(-1, -2).conj())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此，我们会说 `hcore_1` $h_{t \\mu \\nu}^{(1)}$ 是厄米的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 抗磁项"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "抗磁项 $\\hat h {}^{(2)} (\\boldsymbol{\\mathscr{B}}) = \\frac{1}{8} \\big( \\boldsymbol{\\mathscr{B}}^2 \\boldsymbol{r}^2 - (\\boldsymbol{\\mathscr{B}} \\cdot \\boldsymbol{r})^2 \\big)$ 需要一些技巧生成。PySCF 中可以生成张量 `int1e_rr`：\n",
    "\n",
    "$$\n",
    "\\mathtt{int1e\\_rr}_{ts \\mu \\nu} = \\langle \\mu | ts | \\nu \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们定义 `hcore_2` $h_{ts \\mu \\nu}^{(2)}$ 为"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "h_{ts \\mu \\nu}^{(2)} = \\frac{1}{8} \\big( \\delta_{ts} \\langle \\mu | x^2 + y^2 + z^2 | \\nu \\rangle - \\langle \\mu | ts | \\nu \\rangle \\big)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [],
   "source": [
    "with mol.with_common_orig(coord_orig):\n",
    "    int1e_rr = mol.intor(\"int1e_rr\").reshape(3, 3, mol.nao, mol.nao)\n",
    "hcore_2 = 1/8 * (np.einsum(\"ts, uv -> tsuv\", np.eye(3), int1e_rr.diagonal(0, 0, 1).sum(-1)) - int1e_rr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "并且，上述张量具有下述性质：\n",
    "\n",
    "$$\n",
    "h_{ts \\mu \\nu}^{(2)} \\cdot \\mathscr{B}_t \\mathscr{B}_s = \\langle \\mu | \\hat h {}^{(2)} (\\mathscr{B}_t, \\mathscr{B}_s) | \\nu \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Core Hamiltonian 程序实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们最后可以编写外加磁场微扰下的 Core Hamiltonian，以及在此微扰下的分子体系能量。为了加快计算速度，我们会使用为微扰的自洽场密度作为初猜 `dm_guess`。Core Hamiltonian 表达式为\n",
    "\n",
    "$$\n",
    "h_{\\mu \\nu} (\\boldsymbol{\\mathscr{B}}) = h_{\\mu \\nu} (\\mathscr{B}_x, \\mathscr{B}_y, \\mathscr{B}_z)\n",
    "= h_{\\mu \\nu}^{(0)} + \\sum_{t} h_{t \\mu \\nu}^{(1)} \\mathscr{B}_t + \\sum_{ts} h_{ts \\mu \\nu}^{(2)} \\mathscr{B}_t \\mathscr{B}_s\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm_guess = mf.make_rdm1()\n",
    "\n",
    "def hcore_mag_field(dev_xyz):\n",
    "    mf = scf.RHF(mol)\n",
    "    def hcore(mol_):\n",
    "        hcore_total  = np.asarray(scf.rhf.get_hcore(mol_), dtype=np.complex128)\n",
    "        hcore_total += np.einsum(\"tuv, t -> uv\", hcore_1, dev_xyz)\n",
    "        hcore_total += np.einsum(\"tsuv, t, s -> uv\", hcore_2, dev_xyz, dev_xyz)\n",
    "        return hcore_total\n",
    "    mf.get_hcore = hcore\n",
    "    return mf.kernel(dm=dm_guess)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述函数的参数 `t`, `s` 表示坐标方向分量，`dev_t`, `dev_s` 表示外加微扰大小，单位为 a.u.。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "譬如，若在 $x$ 方向的磁场上施加 $\\mathscr{B}_x = 1 \\, \\mathsf{a.u.}$，而 $y$ 方向上施加 $\\mathscr{B}_y = 2 \\, \\mathsf{a.u.}$ (即 $\\boldsymbol{\\mathscr{B}} = (\\mathscr{B}_x, \\mathscr{B}_y, \\mathscr{B}_z) = (1, 2, 0) \\, \\mathsf{a.u.}$)，那么下述程序会给出该自洽场能量 $E_\\mathrm{tot} (\\mathscr{B}_x, \\mathscr{B}_y, \\mathscr{B}_z)$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-57.997213868466154"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hcore_mag_field(0, 1, 1, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数值导数求取磁化率"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们已经有了求取 $E_\\mathrm{tot} (\\mathscr{B}_x, \\mathscr{B}_y, \\mathscr{B}_z)$ 的程序了，接下来就可以进行数值导数计算。数值导数的计算公式可以简单地使用三点差分法；对于被求导量 $x :\\neq y$，有 (当 $h$ 足够小时)\n",
    "\n",
    "$$\n",
    "\\frac{\\partial^2 f}{\\partial x \\partial y} \\simeq \\frac{1}{4 h^2} \\big[ f(x + h, y + h) - f(x - h, y + h) - f(x + h, y - h) + f(x - h, y - h) \\big]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "而对被求导量相同的情形，有\n",
    "\n",
    "$$\n",
    "\\frac{\\partial^2 f}{\\partial x^2} \\simeq \\frac{1}{h^2} \\big[ f(x + h) - 2 f(x) + f(x - h) \\big]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面的程序就依照上述两个公式进行二阶数值导数求取。求导的原点取在 $(\\mathscr{B}_x, \\mathscr{B}_y, \\mathscr{B}_z) = (0, 0, 0)$ 即不受外磁场影响的情形的自洽场能量 `eng_origin`，差分大小为 `interval` $h = 10^{-3} \\, \\mathsf{a.u.}$。需要注意，根据约定俗成，\n",
    "\n",
    "$$\n",
    "\\xi_{ts} = - \\frac{\\partial^2 E_\\mathrm{tot} (\\boldsymbol{\\mathscr{B}})}{\\partial \\mathscr{B}_t \\partial \\mathscr{B}_s}\n",
    "$$\n",
    "\n",
    "因此求取得到的磁化率 `num_polar` $\\xi_{ts}$ 需要乘以 -1。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "eng_origin = hcore_mag_field((0, 0, 0))\n",
    "interval = 1e-3\n",
    "num_polar = np.zeros((3, 3))\n",
    "for t in range(3):\n",
    "    for s in range(3):\n",
    "        if t != s:\n",
    "            dev_xyzs = np.zeros((4, 3))\n",
    "            dev_xyzs[0, t] = dev_xyzs[0, s] = dev_xyzs[1, t] = dev_xyzs[2, s] =  interval\n",
    "            dev_xyzs[3, t] = dev_xyzs[3, s] = dev_xyzs[2, t] = dev_xyzs[1, s] = -interval\n",
    "            num_polar[t, s] = (\n",
    "                + hcore_mag_field(dev_xyzs[0])\n",
    "                - hcore_mag_field(dev_xyzs[1])\n",
    "                - hcore_mag_field(dev_xyzs[2])\n",
    "                + hcore_mag_field(dev_xyzs[3])\n",
    "            ) / (4 * interval**2)\n",
    "        else:\n",
    "            dev_xyzs = np.zeros((2, 3))\n",
    "            dev_xyzs[0, t], dev_xyzs[1, t] = interval, -interval\n",
    "            num_polar[t, t] = (\n",
    "                + hcore_mag_field(dev_xyzs[0])\n",
    "                + hcore_mag_field(dev_xyzs[2])\n",
    "                - eng_origin * 2\n",
    "            ) / (interval ** 2)\n",
    "num_polar *= -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-4.94475,  0.21773, -0.08268],\n",
       "       [ 0.21773, -4.27801,  0.49885],\n",
       "       [-0.08268,  0.49885, -4.15348]])"
      ]
     },
     "execution_count": 178,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_polar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们再与 PySCF 的解析结果作对照："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-4.94475,  0.21773, -0.08268],\n",
       "       [ 0.21773, -4.27801,  0.49885],\n",
       "       [-0.08268,  0.49885, -4.15348]])"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_mag.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[^Atkins-Friedman.Oxford.2010]: Atkins, P. W.; Friedman, R. S. *Molecular Quantum Mechanics*; Oxford University Press, 2010."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
